#include "Multitree.h"
#include <assert.h>
#include <string.h>
#include <gc/gc.h>

Multitree * taylorize(Multitree * tree, Multitree * center, unsigned short int terms, char * variable){
	/*tree must have only one variable (and that variable must be the "variable")*/
	Multitree * x_power = NULL;	/*(x-a)^n*/
	Multitree * XminusA = NULL;	/*x-a*/
	Multitree * result = NULL;		/*whole thing*/
	Multitree * current_derivative = NULL; /*f()'*/
	Number * current_denominator = NULL; /*n!*/
	Multitree * added = NULL;		/*current term*/
	Multitree * temp1 = NULL;		/*Derivative with replaced variable*/
	Multitree * temp2 = NULL;
	Number * numderiv = NULL;		/*Numeric derivative*/
	ListOfChilds * new_child = NULL;
	Number * zero = NULL;
	Number * one = NULL;
	int term_number;
	Number * term_temp = NULL;
	zero = construct_number();
	one = construct_number();
	current_denominator = construct_number();
	term_temp = construct_number();

	inttonum(zero, 0);
	inttonum(one, 1);
	if (tree == NULL || center == NULL || terms == 0 || variable == NULL){
		return NULL;
	}

	/*prepare some helper parts*/

	/*the tree we need to create is + with subtrees in common form: derivative * denominator * (x-a)^power*/
	XminusA = give_me_tree_with_two_incomplete_childs();
	XminusA->symbol.type = OPERATOR;
	XminusA->symbol.oper = PLUS;
	XminusA->childs->inverted = 0;
	XminusA->childs->subtree= treevariable("X"); /*Should not fail - if it does, it means we do not have memory and we can abort*/
	XminusA->childs->next->inverted = 1;
	XminusA->childs->next->subtree = copy_tree(center);

	/*x_power = (x-a)^0)*/

	x_power = give_me_tree_with_two_incomplete_childs();
	x_power->symbol.type = OPERATOR;
	x_power->symbol.oper = POWER;
	x_power->childs->inverted = 0;
	x_power->childs->next->inverted = 0;
	x_power->childs->subtree = XminusA;
	x_power->childs->next->subtree = treefromnum(zero);
	XminusA = NULL;

	/*denominator (0! = 1)*/
	inttonum(current_denominator, 1);

	/*current_derivative*/
	current_derivative = copy_tree(tree);

	
	temp1 = replace_variable_tree(current_derivative, variable, center);
	if(temp1 == NULL){
		return NULL;
	}
	temp2 = canonicalize_tree(temp1);
	if(temp2 == NULL){
		return NULL;
	}
	delete_tree(temp1);
	temp1 = NULL;
	/*try and catch are from C++ version*/
	/*try{*/
	numderiv = numeric(temp2);
	if(numderiv == NULL){
		return NULL;
	}
	/*
	 }catch(...){	//catch everything, destroy structures and throw the exception up
		//we caught either division by zero or unknow variable. In both cases we are not really able to continue
		delete_tree(x_power);
		delete_tree(current_derivative);
		delete_tree(cent_num);
		delete_tree(temp2);
		throw;
	}*/

	delete_tree(temp2);
	temp2 = NULL;
	/*prepare first term*/
	added = give_me_tree_with_two_incomplete_childs();
	added->symbol.type = OPERATOR;
	added->symbol.oper = MULTIPLY;
	/*third part - x_power*/
	added->childs->next->next = give_me_empty_child();
	added->childs->next->next->next = NULL;
	added->childs->next->next->inverted = 0;
	added->childs->next->next->subtree = copy_tree(x_power);
	/*second part - derivative - numeric*/
	added->childs->next->inverted = 0;
	added->childs->next->subtree = treefromnum(numderiv);
	/*firtst part - denominator*/
	added->childs->inverted = 1;
	added->childs->subtree = treefromnum(current_denominator);
	
	/*Added is now the first term we need*/
	if (terms == 1){
		delete_tree(x_power);
		delete_tree(current_derivative);
		return added;
	}
	
	result = give_me_empty_tree();
	result->symbol.type = OPERATOR;
	result->symbol.oper = PLUS;
	result->childs = give_me_empty_child();
	result->childs->next = NULL;
	result->childs->inverted = 0;
	result->childs->subtree = added;

	for(term_number = 1; term_number < terms; ++term_number){
		/*make next derivative*/
		temp1 = diff_tree(current_derivative, variable);
		if(temp1 == NULL){
			return NULL;
		}
		delete_tree(current_derivative);
		current_derivative = canonicalize_tree(temp1);
		if(current_derivative == NULL){
			return NULL;
		}
		delete_tree(temp1);
		temp1=NULL;
		/*increase n in (x-a)^n*/
		numberplus(x_power->childs->next->subtree->symbol.number, x_power->childs->next->subtree->symbol.number, one);

		inttonum(term_temp, term_number);

		assert(numbercmp(x_power->childs->next->subtree->symbol.number, term_temp) == 0);

		/*increase factorial in denominator*/
		numbermul(current_denominator, current_denominator, term_temp);

		/*prepare numeric derivative*/
		temp1 = replace_variable_tree(current_derivative, variable, center);
		if(temp1 == NULL){
			return NULL;
		}
		temp2 = canonicalize_tree(temp1);
		if(temp2 == NULL){
			return NULL;
		}
		delete_tree(temp1);
		temp1 = NULL;

		/*try{*/
		numderiv = numeric(temp2);
		if(numderiv == NULL){
			return NULL;
		}
		/*}catch(...){
			//We delete again everything (including result) and throw exception upwards
			//The code should not get here, because we would get error before (at first call to numeric)
			delete_tree(x_power);
			delete_tree(current_derivative);
			delete_tree(cent_num);
			delete_tree(temp2);
			delete_tree(result);
			throw;
		}
		*/

		delete_tree(temp2);
		temp2 = NULL;


		/*Added (copied code!)*/
		added = give_me_tree_with_two_incomplete_childs();
		added->symbol.type = OPERATOR;
		added->symbol.oper = MULTIPLY;
		/*third part - x_power*/
		added->childs->next->next = give_me_empty_child();
		added->childs->next->next->next = NULL;
		added->childs->next->next->inverted = 0;
		added->childs->next->next->subtree = copy_tree(x_power);
		/*second part - derivative - numeric*/
		added->childs->next->inverted = 0;
		added->childs->next->subtree = treefromnum(numderiv);
		/*firtst part - denominator*/
		added->childs->inverted = 1;
		added->childs->subtree = treefromnum(current_denominator);
		
		/*Now connect it to result*/
		new_child = give_me_empty_child();
		new_child->inverted = 0;
		new_child->subtree = added;
		new_child->next = result->childs;
		result->childs = new_child;
		/*we should probably make result canonical, but I think it is not necessary now*/
	}
	delete_tree(x_power);
	delete_tree(current_derivative);
	temp1 = canonicalize_tree(result);
	result = temp1;
	temp1 = NULL;
	return result;
}

